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In this work we focus on how noise propagates in biochemical reaction networks and affects 
sensitivities of the system. We discover that the stochastic fluctuations can enhance sensitivities in 
one region of the value of control parameters by reducing sensitivities in another region. Based on 
this compensation principle, we designed a concentration detector in which enhanced amplification 
is achieved by an incoherent feedforward reaction network. 



I. INTRODUCTION 

Quantifying biochemical processes at the cellular level 
is becoming increasingly central to modern molecular bi- 
ology [1, [28[ • Much of this attention can be attributed to 
the development of new methodology for measuring bio- 
chemical events in time. In particular, the use of green 
fluorescent protein (GFP) and related fluorophores, high 
sensitivity light microscopy and flow cytometry have pro- 
vided researchers with greater details of the dynamics of 
cellular processes at the level of single cells [H, [28, 38] and 
even single molecules [39j. 

Quantitative laboratory measurements demand theo- 
retical models that are able to describe and interpret the 
data. The traditional approach to modeling cellular pro- 
cesses has been to use deterministic equations with con- 
tinuous dynamics variables: It assumes that concentra- 
tions of metabolites, proteins, etc. vary deterministically 
in a continuous manner. This may be a reasonable as- 
sumption when one is dealing with molecules that occur 
in relatively large numbers. For example, the concentra- 
tion of ATP in Streptococcus lactis is approximately 2.5 
mM [23j which, assuming a cell volume of 10 -15 L, yields 
roughly 1.5 million molecules per cell. In other cases the 
number of molecules can be much lower, for example, 
the number of LacI tetrameric repressor proteins in E. 
coli has been estimated to be of the order of 10 to 50 
molecules per cell (25|. Such small numbers suggest that 
a continuum model may not always be an appropriate 
description. Moreover, given the probabilistic nature of 
chemistry at the molecular level a deterministic approach 
appears to be an inadequate description. In addition re- 
cent experimental measurements have clearly highlighted 
the stochastic nature of biochemical processes in individ- 
ual cells [U . 

Stochastic reaction processes have often been theoret- 
ically investigated using the chemical master equation 
[ill ]. F° r general reaction systems, this equation is a 
challenge to solve analytically because the rates of re- 
actions are often expressed as non-linear functions of 
concentrations (l|. In addition numerical solutions are 
equally impractical because of the huge increase in the 
number of states. Alternative methods such as the Gille- 
spie stochastic simulation algorithm pj]j. can also become 
highly intensive in computation even for reaction systems 



involving several species of molecules. Thus, stochastic 
systems are often modeled with certain approximations 
for analytical and numerical investigations. 

One such approximation is the linear noise approxima- 
tion that has been widely applied to various stochastic 
systems, to estimate variances and co-variances of con- 
centrations [H|]. The mean levels of concentrations are 
predicted to be the same as the conventional determinis- 
tic approach that neglects stochastic fluctuations of con- 
centrations. This approximation however becomes in- 
valid as the numbers of molecules in a system get smaller. 
Each reaction event becomes more distinguishable and 
concentrations fluctuate in greater strengths. The rates 
of reactions, which are dependent on concentrations, be- 
come more affected by the fluctuations. The mean levels 
of the rates can become different from the deterministic 
estimates and this affects the mean levels of concentra- 
tions. The linear noise approximation therefore becomes 
invalid. To correct this discrepancy, different approxi- 
mation schemes have been introduced: mass fluctuation 
kinetics (MFK) [l2j and effective stability analysis (ESA) 
[3lj . These approaches provide more accurate analyses 
at higher noise levels than the linear noise approximation 
approach, by taking into account concentration fluctua- 
tion effects in mean reaction rates. MFK uses the mo- 
ment closure approximation [3o| to describe the evolution 
of the system in terms of the mean and covariance val- 
ues of concentrations in the course of time. It has been 
successfully applied to mass-action type reaction systems 
showing stochastic focusing [27j and noise-induced ge- 
netic oscillators [361 ]. It has, however, some limitations 
on investigating the time evolution of bistable systems 
because mean and covariance values of the concentra- 
tions cannot distinguish bi-stable and mono-stable sys- 
tems [IH . For the study of bistable systems especially in 
stationary states, there exists another approximation ap- 
proach, ESA. By using mode-coupling approximations, it 
successfully describes how concentration fluctuations af- 
fect the bistability. Both MFK and ESA are much less 
intensive in computation than the Gillespie's algorithm 

Here we provide an approximate theoretical analysis 
based on mass fluctuation kinetics to study stationary 
state properties of chemical reaction systems. Our focus 
in this paper is the investigation of how noise propaga- 
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tion [la, |26|, |28|, in different network motifs, affects the 
levels of mean concentrations and mean fluxes. Gomez- 
Uribe et al. [l2j has shown that the difference between 
mean rates of reactions and their deterministic (without 
any noise) estimates can be predicted from the concen- 
tration covariances and the curvatures of the nonlinear 
rate functions (reaction rate equations), as a first order 
contribution [l2j . We have investigated this curvature- 
covariance effect more thoroughly and have found that 
the effect provides simple and clear qualitative illustra- 
tions of stochastic focusing, and leads to qualitative un- 
derstanding on how to design and control chemical reac- 
tion systems to achieve certain noise-responses in system 
behaviors. 

The curvature-covariance effect shows that increased 
sensitivity (stochastic focusing) in one region of the val- 
ues of control parameters can lead to decreased sensi- 
tivity (stochastic defocusing) in other regions. Here the 
sensitivity is a measure of a system response due to a 
source signal change and it is defined as the ratio of the 
percentage change of a response signal (r) to the percent- 
age change of a source signal (s): 
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We have applied this stochastic focusing- defocusing com- 
pensation effect to investi gate an incoherent feedforward 
concentration detector [6|, [2j| • The concentration detec- 
tor shows increased amplification of the concentration de- 
tection. This is due to the stochastic focusing. The sen- 
sitivity of the detection is however not enhanced due to 
the stochastic defocusing. By tuning system parameters, 
we could enhance the amplification up to eight times. 
By applying the curvature-covariance effect, we have fur- 
ther increased the amplification, by modifying upstream 
subnetwork structures. 

We present our analysis in this manuscript as follows. 
In section [II] we will show how both concentration fluctu- 
ations and the curvature of the rate functions affect the 
mean rates of reactions. We also explain the mechanism 
of the stochastic focusing-defocusing compensation. In 
section Hm we will illustrate the fluctuation effects in var- 
ious network motifs, which will also be tested by Monte- 
Carlo simulations by using the Gillespie stochastic sim- 
ulation algorithm. We will investigate negative feedback 
(homeostasis, hyperbolic inhibition) and incoherent feed- 
forward (concentration detection) regulation. 

Our analysis, like other approximation approaches in- 
troduced earlier, is based on the chemical master equa- 
tion and thus the noise is considered intrinsic [35j , which 
means the noise is generated by random chemical reac- 
tions involved in reaction systems and all other noise from 
outside the system is considered negligible. As in MFK 
and ESA, our analysis becomes invalid if the third and 
higher moments of the noise correlation need to be taken 
into account [H, HH ■ Our analysis focuses only on the 
stationary state behaviors of the processes without any 
oscillation, which can be further investigated as an ex- 



tension of this proposed analysis. 



II. RESULTS 

A. Curvature-Covariance Contribution to Mean 
Reaction Rates 

We consider chemical reaction networks where the 
number of particles involved in the reactions are low. Due 
to the small particle numbers, the change in the number 
of particles due to chemical reactions is observed as a 
discrete process [Hj]- In addition, the reactions occur 
due to the random collisions between reactants. Thus, 
the change of molecule numbers is not only discrete but 
also random and is often described by discrete stochas- 
tic processes. Under special conditions (homogeneous 
and statistically independent reaction events), stochas- 
tic processes are fully described by the chemical master 
equation [ll| . This equation describes the time evolution 
of a probability distribution function, which represents 
the probability that one finds the number of molecules 
for each species at a given time. The stochastic reac- 
tion events cause molecule concentrations to fluctuate 
and this determines the rates of the next reaction events. 
These events again cause the concentrations to fluctuate 
and the same argument is applied to the following steps. 
The mean values of the rates can be affected by the con- 
centration fluctuations. This effect will be formulated in 
this section. 

We will characterize the concentration fluctuations by 
using mean values and co- variances of the concentrations: 
E.g., in the fluorescent protein experiments, fluctuations 
in light intensity can be measured in time. Once the 
light intensity fluctuation stabilizes at a constant level 
(in the stationary state), this level indicates the rela- 
tive mean value of the concentration of the light emit- 
ting proteins, and the standard deviation from the con- 
stant level measures the relative strength of the fluctu- 
ations. If the experiments are performed by using both 
green and yellow fluorescent proteins, one can measure 
the intensities of both proteins and quantify how much 
the protein concentrations fluctuate together (fluctuation 
correlation) by measuring their co-variance. 

We consider m species of molecules and n reactions. 
For our purpose we assume that system parameters p 
are non-fluctuating (bold symbols represent matrices and 
vectors), p can be reaction rate constants in mass ac- 
tion rate equations, dissociation constants in Hill equa- 
tions, temperature, pH, etc. The time evolution of mean 
concentrations and concentration co- variances can be de- 
scribed by the following equations [l2|](see Appendix IB) : 
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where Nr is a mo X n reduced stoichiometry matrix 
(l8| . mo is the number of linearly independent rows in 
a stoichiometry matrix, v represents propensity func- 
tions describing the rates of reactions: v = {v%, • • • , u n }. 
We assume that Vi is a function of both {sj} (with 
j = 0, • • • ,mo) and pi, i.e, Vi(s,pi). The angle brack- 
ets denote the average: more precisely, the average of 
data taken at a given time t for repeated independent 
runs of simulations or experiments. J.^ is an element of 
the Jacobian matrix defined as 
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The covariance matrix <x is defined as 



= ((*«-<**»(*i-<*i»)- 



For i — j, the covariance becomes the variance. The 
matrix _D is the diffusion coefficient matrix defined by 
NrAN 1 ^ with a diagonal matrix Ay = ViSij. f2 is the 
system volume. The above equations do not hold ex- 
actly but are valid under certain approximations (see Ap- 
pendix [B] and [Cj . The validity of the above equations, 
however, should be checked on the basis of each different 
model since noise in adjacent systems can be propagated 
into the propensity function [l5|, [2(| [28| . We will show 
several examples of this noise propagation in Example 2. 
We note that Eg. [2] is different from the results of Gomez- 
Uribe, et al. [12 ]: We have neglected all the terms of the 
order of l/f2 2 that have been kept in Eq. 11 of Gomez- 
Uribe, et al. [l2j]. For a detailed discussion readers are 
referred to Appendix [B] and [C] This is consistent within 
the approximation of truncation of third and higher order 
moments. 

The right-hand term in parentheses in Eq. [1] is called 
the mean propensity function and shows that the mean 
rate of reaction, denoted by v, is affected by concentra- 
tion covariances <r and the curvature of the propensity 
function (v), d 2 v/dsidsj. To illustrate this curvature- 
covariance correction, let us consider an example: a sub- 
strate is converted to a product through an enzyme re- 
action and the enzyme-substrate complex turnovers so 
fast that we can assume that the complex is in quasi- 
steady state. Then, the creation rate of the product can 
be described by the Michaelis-Menten rate equation v(s) 



v(s) 



> p. 



s 



We assume that fluctuations in the concentration of 
molecule S, s, are symmetric with respect to its mean 
value. Then, the fluctuation of v becomes asymmetric 
because the positive fluctuations of s cause the fluctua- 
tion of v to become relatively smaller than the negative 
fluctuations of s, as shown in Fig. [IJ due to the curved 
shape of the propensity function v(s). As the curvature 
of the propensity function increases, the probability dis- 
tribution of v becomes more asymmetric and the mean 



value of the propensity function (v(s)) deviates from the 
deterministic prediction v((s)). As the distribution func- 
tion of s gets narrower, the distribution function of v 
also gets narrower and the difference between v((s}) and 
{v(s)) gets smaller. This is why Eq.[T]shows that the first 
order of correction to a deterministic rate is proportional 
to both concentration covariances and the curvature of a 
propensity function. 
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FIG. 1: A mean value of a Michaelis-Menten propensity func- 
tion v(s) depends on its curvature and the variance of sub- 
strate concentration s. The probability distribution function 
of s, P(s), is assumed to be symmetric. The probability dis- 
tribution function of v(s), P(v), becomes asymmetric due to 
the nonlinear propensity function v(s). This makes the dif- 
ference between v((s)) and (v(s)}. 

The non-zero covariances occur in general because the 
concentration fluctuation in one species is propagated 
into that of another species via the connected paths of re- 
action networks and such propagated fluctuations cause 
the propensity function to fluctuate. Such non-local ef- 
fects are incorporated into the covariance term in Eq. [T] 
and we will give an example of this effect in section Iffll 
Example 2. 

We now focus on stationary state properties of Eq. [T] 
and [H In the stationary state, the right hand sides of 
these equations vanish. Eq. [5] leads to a Lyapunov equa- 
tion showing how concentration covariances u are related 
to D (strengths of random driving forces causing the con- 
centrations to fluctuate) and J (strengths of tendencies 
returning to stable fixed points of the case in the absence 
of the random driving forces): 
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The above equation can be considered as a non- 
equilibrium steady-state extension of the equilibrium 
fluctuation-dissipation theorem [2l|, [2(| . By solving the 
above equation for <r, <r can be expressed as a function of 
(s). We denote the solution by u*((s),p). We substitute 
this into Eq. Q] 



N a v({a),p)=0, 



(4) 
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where 



"{(s),p) = v((s),p) 
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2 dsids 



S = (S) 



(5) 

We have denoted the mean propensity function using v. 
The mean propensity function u becomes a true measur- 
able quantity when Eq.|4]is solved for (s) and its solution, 
denoted by s* , is substituted in Eq. [5] 

This shows that the traditional deterministic propen- 
sity function must be corrected in order to yield the 'true' 
measurable mean rate in the stochastic regime by replac- 
ing v(s*,p) by i/(s*,p). This implies that all the theo- 
rems of metabolic control analysis [g, HI, [U Oil , which 
describes how the changes in enzyme activities affect 
metabolite concentration levels, can be applicable to the 
stochastic reaction systems (for more discussion on this, 
we refer to section ITvT). 



B. Stochastic Focusing-Defocusing Compensation 

We investigate the properties of the mean propensity 
function v. Eq. 0J The curvature-covariance effect leads 
to an intuitive understanding of the origin of stochastic 
focusing and furthermore it allows to discover a typi- 
cal feature of the stochastic focusing that the stochas- 
tic fluctuations can enhance sensitivities in one region 
of the value of control parameters by reducing sensitiv- 
ities in another region. We call this stochastic focusing- 
defocusing compensation. This compensation typically 
appears in inhibition regulation and sigmoidal responses. 

Stochastic (de-)focusing is a phenomenon where sys- 
tem sensitivities are enhanced (reduced) due to stochastic 
fluctuations [27|. To understand the stochastic focusing- 
defocusing compensation effect, we consider the following 
reaction system: 



v(s) 



p 



(6) 



where s and p are the concentrations of S and P respec- 
tively and ci and k degradation rate constants of S and 
P respectively, cq is the creation rate of S. We observe 
how a response signal (chosen to be (p)) changes due to 
the perturbation of a source signal ((s)). The sensitivity 
is defined as 



Sensitivity = 



dln(p) _ dln(v(s*)/k) _ dlnu(s*) 
d\n(s) dins* dins* 



where we have used the fact that in the stationary state 
the mean concentration of p is equal to (v(s))/k ~ 
v(s*)/k because the mean creation rate of P balances 
its mean degradation rate. Thus, we investigate how v 
changes due to stochastic fluctuations. 

We consider a sigmoidal response in v(s) given by k\ + 

jj|q^3 with fcj positive constants for i = 1, 2, 3. As shown 
in Fig. [U[b), the curvature of v(s) changes from positive 



sign to negative as s increases from zero. The sign of v— v 
changes from positive to negative due to the curvature- 
covariance effect, and v — v converges to zero as (s) — > 
oo and (s) — > since the curvature vanishes in these 
limits. Therefore, stochastic defocusing (SD) appears in 
between two stochastic focusing (SF) regions as shown in 
Fig.EJb). 

If v(s) represents a hyperbolic-type inhibition of P by 
S: v(s) — fci/(fc 2 + s), the curvature of v is positive for 
all s except s — ► oo. The variance of s vanishes when 
(s) — » 0. Thus, v — v is always positive except that 
(s) = and oo. This means that SD changes to SF as 
(s) increases from zero as shown in Fig. [2ja). 

However, SF does not always come with SD. For ex- 
ample, when v(s) is a Michaelis-Menten type propen- 
sity function, only SF appears without SD as shown in 
Fig.^c). 

We hope that this stochastic focusing- defocusing com- 
pensation effect and the curvature-covariance effect can 
be applied to design and control reaction networks to 
improve system functionality by exploiting intrinsic noise 
[6ll24l|. As an application, in section HTTl we will design in- 
coherent feedforward concentration detectors to improve 
detection amplification by applying our results. 



III. APPLICATIONS 

Example 1: Curvature-covariance Contributions to 
Mean Propensity Functions 

To examine in detail the contribution that the 
curvature-covariance makes to the mean propensity func- 
tion, we will present three examples. In the first example, 
we consider an association reaction [e|: 

i a v=ks 1 s 2 v 
Ol + 02 > A. 

In the stationary state, the association propensity func- 
tion v has its mean value at 

(v) = k(sxs 2 ) = fc(si)(s 2 ) + k{(sx - (si))(s 2 ~ (s 2 ))}- 

The first term on the right hand side is the propensity 
function for deterministic systems. The second term is 
the covariance between s\ and s 2 . Thus, the true mean 
value of the reaction rate can be estimated by taking into 
account the covariance effect. 

As a second example, we consider a Michaelis-Menten- 
type reaction: 
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u(s) 



X, 



(7) 



— Vm 



where vq is a fixed creation rate of S and v(s) = K +s , 
with K m the Michaelis-Menten constant and V max a sat- 
uration rate. If the variable s is non-stochastic, then the 
probability distribution function of s is a delta-function 
centered at (s). Thus, (v(s)) = v((s)). However, if s 
fluctuates stochastically, this equality does not hold any 
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s log(s ; 

(a) Inhibition Regulation 





s log(s ; 

(b) Sigmoidal Response 





s log(s; 
(c) Michaelis-Menten Rate Equation 



FIG. 2: Stochastic focusing-defocusing compensation: Three 
different types of propensity functions for a reaction scheme 
Eq. © show different compensation patterns, s* is a mean 
level estimate of concentration s, which can be varied by 
changing a system parameter such as the creation rate of s. 
Depending on this parameter value, stochastic focusing (SF) 
or stochastic defocusing (SD) appears. Inhibition regulation 
and sigmoidal response of v(s) lead to the compensation [(a) 
and (b)]. Michaelis-Menten type rate equation can result in 
only SF without SD as shown in (c). Solid black lines corre- 
sponds to v(s*) and red lines to u(s*). 



more. The fluctuations in s cause fluctuations in the 
propensity function. Since the propensity is a curved- 
down function in s, the negative direction of fluctuation 
in s will cause the fluctuation in v to be more negative 
than the positive direction of fluctuation in s. Such bi- 
ased fluctuations in v makes (w(s)) smaller than v({s)). 
For small enough fluctuations such differences can be 
shown to be proportional to the curvature of v and also 
to the concentration variance in s: 



1 d 2 v 

2 5? 



(8) 



The above equation is derived from Eq. [5l s* is the 
solution of Eq. fj] and a* the solution of Eq. [3] in the 
stationary state. Such curvature-variance correction to 
the deterministic propensity function is closely related 
to stochastic focusing [27J (see Example 3 and 4). 

Monte-Carlo simulations based on the Gillespie 
stochastic simulation algorithm [Icl ] were performed for 
the above Michaelis-Menten-type reaction. We set V = 
S/ {Km + S) where we rescale time so that the maximum 
rate of V is set to one and varies Vq from to 0.9 for 



given values of Km- An upper case letter S denotes the 
molecule number of species S and the lower case letter s 
its concentration, i.e., s = S/Cl. The propensity function 
V also needs to be divided by the system volume fl to 
become the propensity function v defined in the previous 
section. To avoid too many notations, we take the same 
symbol for both the propensity functions in spite of the 
difference. Km needs to be also divided by fl to become 
the Michaelis-Menten constant. For Km > 1, the vari- 
ance corrections become accurate as shown in Fig. [3] and 
HI For K M = 5 (Fig. Mb)), when v ^ 0.1, S fluctuates 
between and 4 for most of the time in the stationary 
state, where the propensity function is almost linear in S 
(Fig. [3Ja)). Thus, the distribution of S becomes similar 
to the Poisson distribution. This is why the variance cor- 
rection becomes negligible for this range of the value of S. 
For Vq ~ 0.9 the reaction is saturated and the propensity 
function becomes almost linear again. The variance cor- 
rection becomes negligible. For vq values between 0.4 and 
0.7 the variance correction becomes significant but still 
gives reasonable estimates for the mean propensity. For 
the smaller values of Km, the corrections become less ac- 
curate. This is because d 3 v/ds 3 becomes larger and the 
neglected terms in Eq. [8] become significant. Here, we 
have shown the origin of the curvature-covariance effect 
on the mean propensity functions and have also illus- 
trated how accurate the moment closure approximation 
is. The approximation needs to be verified however case 
by case; e.g., the creation rate vq in the above exam- 
ple was assumed to be constant but it can be a function 
of other species concentrations, noise of which, possibly 
significant amounts of noise, can be propagated into vq. 

We are going to investigate reaction systems with non- 
linear propensity functions, e.g., Michaelis-Menten ki- 
netics. This is because various network motifs can be 
directly represented by propensity functions. Such non- 
linear rate equations can be realized in mass-action re- 
action systems, when the time scales of reactions can 
be separable into slow and fast ones and the molecules 
involved in the fast reactions are treated as in station- 
ary states [U, EH HI] • This quasi-steady state approx- 
imation leads to non-linear rate equations with simpli- 
fied network structures. However, this approximation 
results in neglecting all the correlations between the fast 
and slow concentration fluctuations. If the time scales 
are not separable, such correlations become significant 
and the results derived by using the nonlinear propen- 
sity functions can be significantly different from the true 
measured ones. It is however important to investigate 
the simplified network by using the nonlinear propensity 
functions as a first step toward to understand the system 
behavior related to each different network motif. 

To illustrate how reasonable it is to use the nonlin- 
ear propensity functions, we consider a negative feed- 
back system showing homeostasis: A protein species Si 
enhances the phosphorylation of another protein species 
and its dephosphorylated form (S2) inhibits the gen- 
eration of protein Si as shown in Fig. 5(a) We as- 
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FIG. 3: Molecules S are created with a constant rate vo and 
degrades with a Michaelis-Menten rate v(S) = S/(Km + S) 
as shown in Eq. [7] (a) The stationary state degradation rate 
in both the stochastic and deterministic cases are compared 
with numerical simulations for fixed value of Km = 5 for var- 
ious creation rates. The Gillespie stochastic simulation algo- 
rithm is used. (v(S)) and (S) are estimated by time averages 
of a single run of simulation. S* denotes the approximate 
estiamate of (S) given by solving Eq. [3] (b) Probability dis- 
tributions in S are shown for each different value of Vq- 



sume that the phosphorylation-dcphosphorylation cycle 
(v3 and V4) is very fast so that the fluctuation of Si im- 
mediately appears in the concentration of £2 and inhibits 
the generation of Si . We also assume that the number of 
molecules involved in cyclic reactions are large enough 
that the internal noise due to the phosphorylation- 
dephosphorylation is negligible. In this case, the above 
reaction can be further simplified to a negative feedback 
reaction system as shown in Fig. |5(b)| 

This system can show homeostasis due to a strong neg- 
ative feedback. In Fig. [6j the negative feedback becomes 
very strong within a narrow range of Si centered around 
Si = 50; the concentration of Si can be decreased right 
after the degradation is accelerated by increasing P2 but 
the concentration eventually return to the approximate 
value of the original concentration because the creation 
rate is strongly increased. This phenomenon is called 



FIG. 4: Degradation propensity functions in a reaction system 
Eq.[jJfor different values of Km- 



homeostasis. Homestasis has been shown to be re lated 
to suppressing concentration fluctuations (see Fig. |6(b)| 
and |4j). The strong negative feedback corresponds to the 
large difference between unsealed elasticities Q of reac- 
tions vi and V2 ■ Such a large difference means that when 
the concentration fluctuates with respect to the mean 
value of the concentration, the system has a strong ten- 
dency to dampen the fluctuation and return to the mean 
value. When the value of Si is within the range, where 
homeostasis appears (p2 = 1 in Fig. |6(b)[ ) , the probability 
distribution of Si becomes narrow. When p 2 increases to 
3, however, the homeostasis vanishes and the probability 
distribution of Si shows very large fluctuations. 

We compare the above two non-simplified and simpli- 



fied reaction systems, Fig. 5(a) and 5(b) As the cyclic 
reaction gets faster, the distribution function of Si for 
the non-simplified reaction converges to that for the sim- 
plified one as shown in Fig. [7J This shows how reasonable 
the use of non-linear propensity functions is. However, if 
the cyclic reactions are not fast, then the correlation be- 
tween the fast and slow variables Si and S2 (S3) become 
significant and the system behavior becomes very dif- 
ferent from the one with nonlinear propensity functions. 
The variance of Si can increase rather than decrease as 
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p 2 changes from 3 to 1 for the slow reaction case (The 
graph for this result not shown). 

v,=p,S 2 v 2 =P2 S i w 

; — ► S ►A 

7 ■ 1 




(a)A Detailed Reaction Network 

fl( S 1,Pl,P3,P4,P5,Pe) fz=P2Sl 

r ►s, 

I ■ 1 



(b)A Simplified Reaction Network 



FIG. 5: (a) Cycle reactions involving S% and S3 act as a 
negative feedback on Si. The rates of the cycle reactions are 
given as vs = psSiS2/(p5 + S2) and V4 = p4Sa/{p6 + S3) [13]. 
(b) The reaction network shown in (a) becomes simplified 
when the cycle reactions are fast. /1 represents a negative- 
feedback propensity function (for its detail functional form, 
we refer to Tyson et al. HH). A circular (flat) end dotted line 
corresponds to positive (negative) regulation. 



Example 2: Non-local Mean Propensity Functions 

As shown in Example 1, the mean propensity function 
is affected by the noise covariance er. Here we discuss 
the non-local property in the noise covariance. We con- 
sider a negative feedback reaction system as shown in 
Fig. 8(a)| Si produces S 2 , which accelerates the degrada- 
tion of Si . This reaction system can be further simplified 
as in Fig. |8(b)| when the life time of Si is much shorter 
than that of S 2 . The propensity function v\ does not 
depend on k 2 , but its mean value vi becomes dependent 
on k 2 in the stationary state because the concentration 
variance depends on k 2 \ the concentration fluctuations 
are due to the events of both the reactions, so the con- 
centration variance of s depends not only on v\ but also 
on v 2 : 



vi({s), ki, k 2 ) 
ki 



ki(ki + k 2 (s)(K M + (s))) 



K M + (s) 2{K M + (s)) 2 (h + k 2 (K M + <s)) 2 ) ' 

where we have expressed the variance a in terms of (s) 
by using Eq. [3] and [5J 

As a second example, we consider a two-step cascade 
reaction system studied by Paulsson et al. [27| as shown 
in Fig. [51 Signal molecules Si inhibits the production of 
S 2 . The fluctuation in Si is propagated into the propen- 
sity function U3. Thus, the mean propensity function 1/3 
will depend not only on p% and p^ but also on pi and p 2 
due to noise propagation from an upstream reactions vi 
and v 2 . We will study this system in detail in the next 
section. 
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FIG. 6: Negative feedback homeostasis: The number fluc- 
tuation is suppressed by a negative feedback as shown in 
5(b)| The fluctuation strength decreases as the mean 
jer of Si increases by changing P2 from 3 to 1. (See 
Tyson, et al. [33 ] for the relationship between the detailed 
one and the simplified one: pi — 1, ps = 1, pa = 1 and 

s 2 + s 3 = 100.) 
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FIG. 7: The probability distribution of Si for the non- 



simplified reaction Fig. 5(a) and simplified one Fig. 5(b) For 
both the slow and fast cycle reactions, we have used the fol- 
lowing parameters pi — 1, P2 ~ 1, f>5 = 1, P6 = 1 and 
S2 + S3 = 100. For the slow reaction, p 3 — 0.01 and P4, = 0.5 
are used. For the fast reaction, P3 = 1000 and P4 — 50000. 
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(a)A Detailed Reaction Network 



(b)A Simplified Reaction Network 



FIG. 8: A hyperbolic inhibition reaction network: In (a), a 
product of Si accelerates the degradation of itself. In (b), 
this reaction gets simplified into a hyperbolic-type inhibition 
reaction when the reactions involving Si is much faster than 
that of S 2 . fci = P1P2/P4, K M = P2/P4, and k 2 = ps. 
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of Km > 1. However, for lower values of % the correc- 
tion becomes less accurate where the stochastic focusing 
can become significant. This shows the limitation of our 
approximation. We also show the above defined sensitiv- 
ity coefficients for the case of Km = 0.1 in Fig. 
The sensitivity is enhanced for the region of (Si 
while reduced for (S) < 1. Here stochastic focusing in 
one parameter region comes with stochastic de-focusing 
in another region. This also can be understood by us- 
ing the curvature-covariance correction effect in the mean 
propensity function 1/3. At (S) = 0, there is no variance 
and thus V3 — V3. At (S) = 00, the hyperbolic curve 
becomes linear and its curvature vanishes, resulting in 
z/ 3 = v$. As (S) decreases from 00, V3 increases faster 
than i>3 but as (S) becomes closer to zero, V3 increases 
slower than V3 (see Fig. 11(b)). Thus, stochastic focus- 
ing and defocusing appears one after the other. This 
implies that one can achieve higher sensitivities of one 
region by sacrificing sensitivities of another region. Wc 
use this stochastic focusing-defocusing compensation ef- 
fect to enhance the amplification of concentration detec- 
tion designed by using incoherent feedforward networks, 
in the next example. 



FIG. 9: Stochastic focusing and defocusing by using a two- 
step cascade reaction system, where Si inhibits the creation 
ofS 2 . 



Example 3: Stochastic Focusing-Defocusing 
Compensation 

The curvature-covariance effect can explain stochastic 
focusing, which describes the phenomenon that system 
sensitivities are enhanced due to stochastic fluctuations 

We consider again the example of the cascade reaction 
as shown in Fig. [5J In Paulsson et al. .27], stochastic 
focusing was studied by perturbing the parameter p\ and 
examining how sensitive the mean concentration of S2 is 
to the change in the mean concentration of S\. We can 
quantify such sensitivity in terms of concentration control 
coefficients: 



where C£ represents the percentage change of x due the 
percentage change of a parameter p from one stationary 
state to another. 

The mean propensity V3 becomes larger than the de- 
terministic rate V3 because the curvature of V3 is posi- 
tive. Since 1*4 is a linear function to S2 (no curvature), 
j/4((s)) is the same as Ui((s)). In the stationary state, 
the mean degradation rate V4 becomes balanced with V3. 
Thus, the mean level of S2 will be enhanced as shown in 
Fig. 1101 The mean propensity function v estimates the 
true mean rate of reaction quite accurately for the value 



The stochastic focusing we discuss has minor concep- 
tual differences from the one discussed in Paulsson et 
al. [23. There stochastic focusing is the only effect of 
a rapidly fluctuating Si. However, here the stochastic 
focusing is independent of how fast Si fluctuates. This 
difference is due to the fact that the focus in Paulsson 
et al [27j was on what is the most probable state of S2, 
and we focus on the mean value of the Si. To under- 
stand this difference, we need to understand the dynam- 
ics of S2 that is correlated with the fluctuation in Si. 
When all reaction rates v\, V2, V3, and V4 are in the same 
order of magnitude, Si and S2 fluctuate on the similar 
time scales. If Si hits zero the inhibition of S2 is re- 
moved and thus the number of S2 can rapidly increase to 
a very large number. When Si increases to 1, however, 
the inhibition acting on S2 appears and the number of 
S2 rapidly decreases. Therefore, the time series profile 
of S2 shows a flat lower bound at zero with many large 
sharp spikes. This time series profile changes as the cre- 
ation and degradation of Si get faster. If the parameter 
values of p\ and P2 increase such that Si fluctuates much 
faster than S2, S2 sees the averaged behavior of Si and 
is unlikely to hit zero. Thus, the strong/weak inhibition 
by Si is averaged out. Due to this averaging, the sen- 
sitivity enhancement, i.e., stochastic focusing, manifests 
itself. This is why the stochastic focusing was claimed to 
be the effect of a rapidly fluctuating Si in this example 
[27j . However if one can observe the time series profile 
for a very long time to get good statistics of the spike 
heights, the mean value of S2 (the mean propensity V4) 
becomes actually independent from how fast V\ and V2 
are, if the ratio P2/P1 is presumed to be kept constant. 
This is why the stochastic focusing defined here becomes 
time-scale independent. 



9 




10"" 10"' 10 

S 1 , <S 1 : 
(c)K M = 0.01 



10' 



10" 



FIG. 10: The stationary state degradation rate of S2 vs. 
stationary state mean number of Si in a cascade reaction 
system as shown in Fig. [9] p 3 = 10000, Km = 1,0.1,0.01, 
and P1/P2 are varied from 1 to 100. 



Example 4: Concentration Detection from 
Incoherent Feedforward Networks 



In this section we investigate how the stochastic 
focusing-defocusing compensation effect can be directly 
related to the sensitivity change of an incoherent feed- 
forward network acting as a concentration detector [f|. 
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FIG. 11: (a) Sensitivity control coemciet (Eq. [9| vs. mean 
molecule number of 5*i for a cascade reaction system as shown 
in Fig. [9]for Km = 0.1. P1/P2 are varied from 1 to 100. (b) 
Stochastic focusing-defocusing compensation: Mean values of 
a propensity function V3 becomes larger than the estimate in 
the deterministic case due to the curvature- variance effect. 



The stochastic focusing-defocusing compensation effect 
will be shown to explain why stochastic fluctuation can 
amplify concentration detection while the sensitivity of 
detection is not enhanced. 

We explain first why incoherent feedforward networks 
can result in concentration detection. As an example, 
consider a species (S2) is regulated by two different path- 
ways: either directly by Xq or indirectly via Si as shown 
in Fig. [T3J The direct control acts as an activator for 
the production of S2 while the indirect control acts as an 
inhibitor (see Fig. 112]). Thus, the feedforward is called in- 
coherent. When the concentration of Xq is zero, S2 is not 
created. As Xq increases, S2 increases together but when 
Xq becomes larger than a thresh-hold point it begins to 
decrease and eventually is dominated by SVs inhibition 
(Fig. flU]) . Thus, one can detect a specific range of the 
concentration of Xq by monitoring the concentration of 
S 2 . 

As shown in Fig. [T3] concentration detection is am- 
plified compared to the deterministic case, but sensi- 
tivities are not enhanced. It is due to the stochastic 
focusing-defocusing compensation effect. The networks 
represented in both Fig. [5] and 12(a) look identical ex- 
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X and Si are not changed (see Fig. [12] (c)). Although 
the variance of Si increases due to the noise propagated 
from Xo, the amplification is reduced. This is because 
we have not taken into account the variance effect of Xo, 
which contributes to the amplification negatively. From 
Eq. [5] the curvature-covariance correction term is given 
by 



1 d 2 v 
2d^ 1 



+ 



d 2 v 
dxodsi 



ld 2 v 

2d^ c 



The first (third) term is negative (positive) because 
the curvature is negative (positive) with respect to the 
change of £0 (si)- The second term vanishes because the 
covariance between Xo and Si vanishes. The vanishing 
covariance happens to be true only in the stationary state 
0, H3, H3| ■ Thus the amplification becomes smaller than 
the original case. 



IV. DISCUSSION AND CONCLUSION 



FIG. 12: Incoherent feedforward reaction systems: S2 is cre- 
ated from Xq under the inhibition control of Si that also be 
created from Xo. Xo is a boundary species of which the con- 
centration is not allowed to fluctuate in (a) and (b). In (c), Xo 
is allowed to fluctuate, (a) The fluctuation in Si follows the 
Poisson distribution, (b) The fluctuation in Si is wider than 
the case (a) with the mean concentration of Si remaining the 
same, (c) Xo is allowed to fluctuate. 



cept that Si and S2 are created from the common source 
in Fig. 12(a) Thus, stochastic focusing appears for large 
values of (Si), which means for large Xo, and stochas- 
tic defocusing for the smaller values. In the case of 
Fig. 12(a) the detection can be enhanced by approxi- 
mately 30% compared with the deterministic case. For 
a different parameter set as shown in Fig. p(b)| the de- 
tection can be enhanced by more than eight times. This 
is because this parameter set gives a similar condition to 
the case of the cascade reaction network for Km = 0.01 



as shown in Fig. 10(c) where stochastic focusing becomes 
significant. 

To enhance the amplification, we exploit the curvature- 
covariance effect. Since the curvature of W3 with respect 
to Si is positive, stochastic focusing gets stronger with 
the increase of the variance of S\. To increase the vari- 
ance, we replace the upstream reaction network of the 
creation and degradation reactions of Si as in Fig. fTJ] 
(b): Xo creates two Si molecules with the same reaction 
rate and degrades two times faster. Thus, the mean val- 
ues of the concentration of Si docs not change but its 
fluctuation is shown to increase |4(. This further ampli- 
fies the detection as Fig. [TBI 

We modify the original upstream network to increase 
the noise fluctuation of Si in a different way. We allow 
Xo to fluctuate, while both the mean concentrations of 



We have analyzed the stationary state properties of 
stochastic reaction systems based on the mass fluctua- 
tion kinetics focusing on how intrinsic noise prop- 
agation affects system sensitivities. We have considered 
nonlinear propensity functions such as Michaelis-Menten 
rate equations as a first step toward to understanding the 
relationship between network topologies and sensitivity 
changes. 

We have investigated how mean levels of concentra- 
tions can be estimated by using the correction to the 
deterministic prediction by taking into account both cur- 
vature effect of propensity function and concentration co- 
variances. The curvature-covariance correction has been 
applied to predict stochastic focusing, concentration de- 
tection, and bistability (readers are referred to Appendix 
ID]) qualitatively. Our analysis shows that stochastic fo- 
cusing comes with stochastic de-focusing typically in the 
systems showing sigmoidal responses in the propensity 
functions and hyperbolic-type (zero order unltrasensitiv- 
ity) negative feedback. As an application of the compen- 
sation effect, we have investigated incoherent feedforward 
concentration detectors. The exploitation of the stochas- 
tic focusing can leads to the amplification of the concen- 
tration detection. The detection sensitivity is however 
not enhanced due to the stochastic de-focusing. We fur- 
ther analyzed how the amplification is affected by addi- 
tional intrinsic noise. We have provided two cases that 
the amplification is enhanced and decreased and have 
explained the behavior by using our curvature-covariance 
correction. The upstream sub-network structure is shown 
to be significantly important to control the amplification. 
We have also analyzed a bistable reaction system having 
a positive-feedback reaction network. We have explained 
why intrinsic noise can reduce the bistability by using 
our analysis (see Appendix [D)) . 

As a further application of our analysis, we have shown 
that the stochastic version of MCA theorems exists: these 
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(b)Significant Amplifications Due to Stochastic Effects 

FIG. 13: A Concentration Detector: mean concentration be- 
comes very large only within a narrow region of the value of 
Xq. Three different cases shown in Fig. [12] are compared. For 
the case (a), the estimates based on our analysis [Curvature- 
Covariance: Case(a)] is compared with its determinstic case 
and its stochastic simulations [Gillespie: Case (a)]. Parame- 
ters: po = 1, pi = 1, P2 = 1, P3 = 0.01, Pi = 1, ps = 0.001 
for subfigure (a). Significant amplification can be achieved 
more than 8 times of the deterministic case as shown in sub- 
figure (b). Parameters: pi — 1, P2 — 1, Pz = 0.01, P4 = 100, 
Ps = 0.001 



theorems can be derived by replacing v in MCA theorems 
to v. We name this extension stochastic control analy- 
sis (SCA). However, SCA has one drawback: elasticity 
(of a mean propensity function v) does not become a lo- 
cal sensitivity measure any more. This is because noise 
propagation can affect the mean propensity function. In 
MCA, the global sensitivities (control coefficients) are ex- 
pressed in terms of the local sensitivities (elasticities) and 
the system-wide response can be described by the com- 
bination of the local response. However, in SCA such 



a description is not possible due to the non-local prop- 
erties of the elasticities. To resolve this issue, the noise 
propagation needs to be described by modular structures, 
i.e., local transfer functions 0, [H, [37| • We have been 
investigating the modular structure of the noise propaga- 
tions as our current research to convert SCA into a useful 
form. For a quantitative analysis, SCA needs to be ap- 
plied to mass-action reaction systems without assuming 
quasi-steady state approximation [l2j ■ 

Our analysis is focused on the stationary state static 
responses. However, it can be extended for the study of 
dynamic responses by taking Fourier transformations of 
Eqs.[T]and[2]just like dynamic res pon se studied in classic 
metabolic control analysis in [lrl |30| . One of the appli- 
cations of the dynamic response analysis is feedforward 
networks acting as frequency filters. 

We hope that the curvature-covariance effect and the 
stochastic focusing-defocusing compensation effect can 
be used for intuitive understanding of how stochastic in- 
trinsic noise affects the system functions. This under- 
standing will help to systematically design and improve 
a synthetic genetic circuit 0, [f| , especially how to control 
the mean levels of concentrations by modifying subnet- 
works in stationary states. 
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APPENDIX A: 



THE CHEMICAL MASTER 
EQUATION 



Chemical reactions occur in a random fashion due to 
collisions between reactants. We assume that each reac- 
tion event arises homogeneously (in a uniform fashion in 
position space) and also independently to another event. 
Such reaction systems are often modeled by stochastic 
processes. The processes are fully described by the time 
evolution of probability distribution of reactant number 
for each species at a given time, mathematically formu- 
lated by the chemical master equation [111 ] . This equation 
describes the time evolution of a probability distribution 
function, which represents the probability that one finds 
the number of molecules for each species at a given time. 

We consider m species of molecules involved in n reac- 
tions: 



Vi 



mi Si 



+ m m S. 



where the molecule numbers are denoted by {Si} with 
i = 1 , ■ ■ • , in and the rate of reaction by Vi . We also 
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assume that Vi can be controlled by changing a parame- 
ter pi. The number change in species i by a single event 
of the above reaction / is described by a reduced stoi- 



chiometry matrix: N R . 



= ml 



. The numbers evolve 



stochastically and their evolutions can be described by 
the chemical master equation: 



dP(S,t) 

at 



E [p({Si - N R J,t)V 3 ({S t - N Rij , Pj }) 

3=1 

-P(S,t)^(S jft -)l. (Al) 



The first term in the right hand side corresponds to 
the probability increase due to the state change: {Si — 
N Rij } — > S through events of reaction j occurred at time 
t. The second to the probability decrease due to the state 
change: S -> {& + ATr, }. 



APPENDIX B: DERIVATIONS OF EQS. Q] AND [2] 

We switch the representation of states from numbers 
{S} to concentrations of molecules {s = S/f2} with 
a system volume, since this concentration representation 
has a direct correspondence to deterministic macroscopic 
kinetics. 

The mean concentration of a species i is given as 
(si) = [IfJLi J™= dsj]siP(s,t) where P(s,t) is the 
probability distribution function of molecule concentra- 
tions, s — {si, S2, • • • , s m }. The evolution of the mean 
concentration of a species i is shown later in this section 
to be governed by the following equation: 



d(f> 
dt 



N R (v(s,p)}, 



(Bl) 



where TV^ is a reduced stoichiometry matrix and v rep- 
resents propensity functions: v = {v±,--- ,v n } with 
Vi = Vi/Cl. We note that Vi is a function of {sj} with 
j = 0, • • ■ , mo and p,, where mo is the number of linearly 
independent rows in a stoichiometry matrix. 

A concentration covariance between two species (i and 
j) is defined as 



The correlations between different molecular species i 
and j and the variance of a species i are quantified by er^ 
and <7,j, respectively. The correlation measures the sta- 
tistical independence between two random fluctuations in 
each different species. If the correlation vanishes, their 
concentration fluctuations are independent statistically. 
E.g., if cry = with i =/= j, the fluctuation with respect to 
the mean value of Sj statistically is not related to the fluc- 
tuation of Sj , although the mean value of Si may change 



depending on the mean value of Sj through Eq. IBll The 
variance of Sj , uu , quantifies its fluctuation strength. The 
evolution of the covariance matrix is shown later in this 
section to be described by the following equation: 

^ = ((N R v)(s-(s)) + (s-(s)f (N R vf+°), (B2) 
where the diffusion coefficient matrix D is defined by 



N R AN R with a diagonal matrix Aj 



It is almost impossible to solve equations (|B1|) and 
(1B2[) unless the propensity function v(s,p) is linear with 
s. E.g., for a nonlinear function v(s, e) = s 2 , equation lB2l 
cannot be solved unless the third moment of s — (s), 
((s — (s)) 3 ), is known already. To find the value of the 
third moment, the fourth moment needs to be known 
too. The same argument applies to all the higher mo- 
ments. Thus, we need to truncate the moment series to 
solve these equations [35j . In the following paragraphs, 
we assume that the third and higher order moments are 
neglected, and then we derive Eqs. Q] and [H which de- 
scribe the approximate evolutions of the mean values and 
covariances of the concentrations (l2] |. 

First, equations lBll and lB2l are derived from the master 
equation, Eg. IA1I Equation IB II is derived as follows. 



dt 



5 j=l i=l 
n m 

EEh-^)(n 



VjP 



E 



VjP-SkVjP 



S j=l i=l 
n m 

L v » 

E -N R 

s i=i i=i 



EE^.(n^ 

S j=l 



n m 

EE[(n- 



SkVjP-SkVjP 



S 3=1 

where E~ NRi i is a raising/lowering operator: 
E- NR ^f{S) = f({Sij - N Rtj }). The first term in 
the right hand side vanishes. We now switch this num- 
ber representation to the concentration representation 
by replacing S —> ils and V — > fto. Then, equation IBll 
is derived after the cancellation of f2 in both hand sides 
of the above equation. 

Equation IB2I is derived as follows. First we define, a 
number covariance matrix S: 

The time evolution of this covariance matrix is given by: 
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n m 



dt 



- 1 



S 3 = 1 <=1 

n m 

ee [(n^^)^ - ( s ")+ n *«k s i ( s i)+ n *») ^ (^))(^ - &)) 

S 3 = 1 »=1 
n 

E ((& - < s *> + jwcs - <$) + N *n) v i - - ( s k))(Si - (Sim) 

i=i 

n 



3=1 



(S - (S)Y (N rV) + (N R V)(S - {S)) + N R A N R 

I 



where AL = ViSij. By switching the number representa- 
tion to the concentration representation (S = Cl 2 cr), we 
derive Eq. [R2l 

Equations [1] and [5] are derived from Eqs. lBll and lB2l bv 
using the Tayler expansion of the propensity function v 
with respect to (s) and by neglecting the third and higher 
moments. We will discuss this approximation further in 
Appendix [Cj Then, equation IB 1 1 can be expressed as 



dt 



N R (v((s),p))+Y, 



dv 



S=(S) 



{Si - (Si)) 



d 2 v 



1 m 

- T 

1,3 = 1 ■> 



S = (S) 



(si - (si))(sj - (sj)) 



The second term in the right hand side vanishes because 
(si — (si)) = (si) — (si) = 0. Equation [1] is derived. 

In the similar way, Eq. [5] also can be derived from 
Eq. IB21 The propensity function v is expanded by using 
the Taylor expansion. Then, the first term in the right 
hand side of Eq. IB2I becomes, after neglecting the third 
and higher moments: 

((N R v)(s (*))) = (N R v((s))((s (*))) + J<x, 

where J is a Jacobian matrix defined as Jij = 
Sfe N Rik dvk/dsj\s=(s) ■ The first term in the above van- 
ishes. We take the same procedure for the second term in 
the right hand side of Eq. [52 The third term of Eq|B2l 
has an extra factor of 1/fi and we neglect the terms of 
the order of 1 /CI and the higher in the Taylor expansion 
of the propensity function which comprises the diagonal 
elements of A: 



(Ay) = Sijivi) = Sij Vi((a)) + 



dvj 
~ds~ 



S = (S) 



The second term vanishes. Therefore, equation [5] is de- 
rived. 



APPENDIX C: MOMENT CLOSURE 
APPROXIMATIONS 



In this section, we will explain in detail the approxi- 
mation we have taken to derive Eqs. [1] and [2 The first 
assumption that we have taken is that intrinsic noise is 
not strong enough that the propensity function v can be 
expanded in Taylor series. The Taylor expansion of v(s) 
with respect to (s) should converge in the region that the 
value of s is sampled with high probability. The second 
assumption is that the third and higher moments of s—(s) 
are negligible. We will explain the second assumption in 
detail in the following paragraphs. 

We consider a simple reaction system: 



S 



where the creation rate is a constant. The mean propen- 
sity of v is Taylor-expanded and the second and third 
terms of the Taylor expansion are compared with each 
other. The third is required to be much smaller than the 
second for the approximation to be valid: 



l d 2 v 
2cV 



1 d 3 v 



S =< S ) 



((s-(s)) 3 ). 



If the mean concentration of s is fixed while a system 
volume CI increases, v((s}) is invariant under the change 
of the system volume. Thus, the inequality in the above 
is determined by ((s — (s)) 2 } and ((s — (s)) 3 ). Then, why 
can it be reasonable that ((s — (s)) 3 ) is much smaller 
than ((s — (s)) 2 ) for the large system volume? We as- 
sume that a system is near the stationary state. Oth- 
erwise, the transient process to the stationary state can 
be arbitrary because the process depends on the initial 
condition of the process. If v is almost a linear function 
of s, the distribution of s becomes similar to the Poisson 
distribution. Then, the higher the moments the smaller 
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in magnitude; 



<*» 2 > 
<s» 4 > 



(5)/fi = 0(l), 
((5-(5)) 2 }/0 2 ^ 

<(S-(S» 3 >/n 3 = 
«S) + 3(S) 2 )/0 4 



<S}/0 2 = 0(1/0), 

(S)/n 3 = o{i/n 2 ), 
= o(i/n 2 ). 



If is large enough that the above criteria satisfied, then 
the moment closure approximation is justified for the pro- 
cess similar to the Poisson process. If v is highly nonlin- 
ear in s, the above scaling relationship does not hold any 
longer. If we assume the fluctuation of s is small enough 
that v can be linearized with respect to the mean con- 
centration of s, the distribution of s becomes similar to 
the Gaussian distribution. For the Gaussian distribu- 
tion, one can show the relationship between the second 
and fourth moments: ((s - (s)) 4 } = 3((s - (s)) 2 ) 2 . This 
means that the fourth moment becomes much smaller 
than the second moment if the second moment is very 
small, and the moment closure approximation can be- 
come valid. 

There is another way of power counting. Consider the 
case that the system volume fl is fixed while the molecule 
number increases and the stochastic variable S follows 
the Poisson process. If the degree of the moment of 
(S — (S))/(S) gets higher, the magnitude of the moment 
becomes lower; 

((5-(5}) 2 )/(5) 2 = 1/(5), 
((S-(S)) 3 )/(Sf = l/(5) 2 , 
<(S-(S» 4 )/<S) 4 = l/(5) 3 + 3/(5) 2 . 

(S — (S))/{S) means the percentage change in the 
molecule number from the mean value. Thus it is more 
directly applicable to biological experiments and also 
simulations based on the Gillespie algorithm. For exam- 
ple, consider in GFP experiments one wants to study the 
protein number statistics. Then, the percentage change 
in GFP's can be estimated by measuring the light inten- 
sity fluctuation. In Monte-Carlo simulations using the 
Gillespie algorithm, the simulation does not have volume 
as its parameter; the volume information is rather hidden 
in the rate of reaction or in the reaction time scale. 

For the mathematical ease, we have used the former 
way of power counting to derive the Eq.[T]and[2](see Ap- 
pendix |B|) . To compare with Monte-Carlo simulations us- 
ing the Gillespie algorithm, we switch back to the number 
representations to compute the means and covariances in 
molecule numbers. Thus, the latter way of power count- 
ing applies for Monte-Carlo simulation results. 



APPENDIX D: BISTABILITY FROM A 
POSITIVE FEEDBACK NETWORK 

In this section we investigate how a bistable system is 
affected by concentration fluctuations. We show that the 



concentration fluctuations can reduce the bistability due 
to the curvature-covariance effect. 

The mass fluctuation kinetics (MFK) has been applied 
to a bistable system in Gomez- Uribe et al. [H| • MFK de- 
scribes the time evolution of the system variables (mean 
and covariance values of concentrations) deterministi- 
cally and it cannot describe stochastic switching: jump- 
ing from one local stable state to another. However, the 
existence of the bistability can be determined by solv- 
ing Eq. 21 although in a very qualitative level. When 
concentrations fluctuate around only one of the locally 
stable regions for most of the time and rarely jumps to 
the other stable one, Eq. [4] can provide a correction due 
to the local concentration fluctuation (however, without 
taking into account the fluctuations between the two lo- 
cal stable regions). 

We consider a positive feedback reaction network as 
shown in Fig. [TJ] In the deterministic case, the station- 
ary state value of S can take two different stable levels. 
One corresponds to a low number state and the other 
to a high number state. Now, we take into account 
the concentration fluctuation and we examine how the 
fluctuation affects bistability. In the stochastic frame- 
work, bistability means the existence of double peaks in 
the concentration probability distribution function (see 
Fig. 115( a)). There is always a finite probability that the 
number S meanders back and forth between the low and 
the high number peaked regions (stochastic switching) . If 
the two peaks are well separated, such switching however 
rarely happens. In this case, the probability distribution 
function of S can be considered a superposition of two of 
each mono-stable distribution functions centered at each 
peak. Each distribution gives an estimate of the mean 
concentration level of one of the bistable states. This is 
why our analysis based on the curvature-covariance ef- 
fect can estimate when the system becomes bistable, al- 
though in a qualitative level since the stochastic switch- 
ing is neglected. In this example as shown in Fig. [TJ] the 
curvature-variance correction reduces the bistability; v\ 
larger than v\ for the positive curvature region of v% and 
also makes v\ smaller than v% for the negative curvature 
region. Thus, bistability appears at the larger value of P2 
(see Fig. [T4|) and disappears at its smaller value as shown 
in Fig. U5f bV However, there exist singular points in v\ 
because a* (the solution of Eq. [3]) is undetermined when 
the local slope of v\ and V2 are the same (see Fig. 15(c) |. 



For our analysis to be valid, the mean concentration es- 
timates (the solution of Eq. 0J) needs to be far away from 
the singular points. Such singular behavior explains why 
more than one unstable stationary level of S is predicted 
as shown in Fig. [TSTb) . The prediction of the unstable 
stationary levels should be neglected because the approx- 
imation used in deriving Eq. [3] becomes invalid. 

We have shown that intrinsic noise can reduce bista- 
bility resulting from a positive feedback by considering 
the curvature-covariance effect. Although this analysis 
based on the curvature-covariance effect needs to be per- 
formed with a significant care when applied to a bistable 
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system, a qualitative and intuitive level of understanding 
will be helpful to predict the stochastic noise effect on 
the system. 

s 3 

v i = ^r— +P 2 v 2 =p 3 S 

1 " 

■ ■ 



FIG. 14: A positive feedback reaction system: S itself accel- 
erates the creation of itself and this is represented by the Hill 
function. 
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FIG. 15: Bistability of a positive feedback system shown in 
Fig. 1141 The probability distribution function for S shows 
double peaks depending on the value of P2 in (a) . Stationary 
state values of S* are estimated from each individual deter- 
ministic and stochastic approaches in (b), where solid lines 
correspond to stable stationary stationary state concentra- 
tions and dotted lines to unstable ones. In (c), the propensity 
function, v\, diverges. Two purple lines indicate the stable 
concentrations. (Parameters: p\,pz = 100000, 0.0115) 



